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NOMENCLATURE 


c 

P 

c 

PO 

F 

G 


H 


H 


rb 

u 

v 


X 

4 

5 
n 
♦ 
♦ 
0 


static  pressure  coefficient 
total  head 

reduced  vorticity  —  defined  by  Eq.  (3) 

perturbation  stream  function  —  defined  by  Eq.  (2) 

derivative  of  G  with  respect  to  n  —  see  Eq.  (29) 

3G 

spline  derivative  approximation  of  -5— 

an 

3H 

spline  derivative  approximation  of 
number  of  nodes  in  the  ^-direction 
number  of  £  nodes  on  the  body 
number  of  nodes  in  the  n-direction 
radial  coordinate 
body  radius 

velocity  component  in  x-direction 
velocity  component  in  r-direction 
axial  coordinate 
vorticity  magnitude 

transformed  coordinate  along  body  and  centerline 
transformed  coordinate  away  from  body  and  centerline 
stream  function 

angle  of  tangent  to  body  surface 
meridian  angle 


All  other  quantities  are  defined  in  the  text 
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I.  Introduction 

The  traditional  approach  to  solving  the  axisymmetric  strong  interaction 
problem  is  the  displacement  body  method  of  which  numerous  examples  have 
appeared  in  the  literature  [1-3].  In  this  procedure  the  boundary-layer  dis¬ 
placement  thickness  is  added  to  the  body  to  obtain  the  strong  interaction 
pressure  distribution.  Since  the  displacement  thickness  depends  on  knowing 
the  pressure  distribution  iteration  is  required  to  find  the  solution.  In 
the  direct  mode  (pressure  distribution  given)  the  boundary  layer  solution 
inevitably  produces  a  displacement  thickness  which  has  numerical  noise  in 
the  strong  interaction  region.  This  "noise"  is  then  amplified  into  pressure 
wiggles  in  the  potential  flow  solution.  Thus  to  obtain  convergence  of  the 
solution,  numerical  smoothing  of  the  pressure  distribution  is  required  between 
each  iteration.  This  process  usually  requires  intervention  by  the  user  and 
is  difficult  to  automate.  A  variation  of  the  displacement  body  method  is  to 
solve  the  boundary-layer  portion  in  the  inverse  mode  which  results  in  a 
numerical  smoothing  effect.  An  example  of  this  approach  is  the  work  by 
Carter  and  Wornom  [4] .  Whether  this  approach  can  be  truly  automated  remains 
to  be  seen. 

An  alternate  method  of  treating  the  strong  interaction  problem  is  through 
the  frozen  vorticity  approximation.  Geller  [5]  has  used  this  approach  to 
solve  for  the  velocity  profile  development  in  the  tail  region  of  a  body  of 
revolution.  In  his  procedure  a  local  body  oriented  cartesian  coordinate 
system  is  used  and  3v/3x  is  neglected  in  the  vorticity.  As  a  result  the 
velocity  profile  at  a  given  body  station  is  obtained  by  solving  an  ordinary 
differential  equation.  He  matches  the  inviscid  rotational  profile  with  a 
one-seventh  power  law  turbulent  boundary  layer  representation  valid  in  the 
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law-of-the-wall  region.  Geller  states  that  the  reason  for  the  success  of 
his  method  is  that  the  dominant  mechanism  in  such  a  flow  is  the  radial 
contraction  of  the  material  vortex  filaments  as  they  approach  the  tail. 

In  other  words,  the  pressure  distribution  in  the  strong  interaction  region 
is  primarily  controlled  by  the  vorticity  in  the  outer  portion  of  the 
boundary  layer  where  usually  the  diffusion  length  is  much  greater  than  the 
strong  interaction  length. 

This  report  presents  a  frozen  vorticity  treatment  of  the  axisymmetric 
strong  interaction  problem  but  differs  from  that  of  Geller  in  two  respects  — 
first,  a  Poisson  equation  is  solved  for  the  rotational  flow  field  and  second, 
no  matching  is  performed  close  to  the  wall.  Instead,  a  slip  velocity  is 
allowed  at  the  surface  which  necessarily  must  exist  because  the  viscous  terms 
have  been  dropped.  Also,  the  frozen  vorticity  solution  is  used  only  to 
obtain  the  strong  interaction  pressure  distribution. 

The  main  reason  for  choosing  the  frozen  vorticity  approach  is  to  achieve 
a  high  degree  of  automation  and  thus  to  produce  an  easy  to  use  code.  This 
objective  is  achieved  by  considerable  commonality  in  the  computational  grids 
used  in  the  potential  and  frozen  vorticity  flow  fields.  In  addition,  the 
boundary  layer  solution  uses  the  same  streamwise  step  size  distribution  as 
the  potential  and  frozen  vorticity  solutions.  This  commonality  together 
with  computer  generated  data  files  allows  the  method  to  be  as  user 
independent  as  possible. 
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II.  Analysts 


Equation  of  Motion 


In  cylindrical  coordinates,  (x,r,0)  where  3/39  -  0  (axial  symmetry) 

inviscid  potential  or  rotational  flow  is  governed  by  the  equation  [6] 

G  +  G  -  -  G  =  -  r2F(«)  ,  (1) 

xx  rr  r  r 

where  a  lower  case  subscript  denotes  partial  differentiation,  G  is  the 
perturbation  stream  function  defined  by 

G  =  *  -  I  r2  ,  (2) 

and  F  is  the  reduced  vorticity  related  to  the  vorticity  C  by 

C  =  rF(«J>)  .  (3) 

The  quantity  F  is  uniquely  determined  by  the  upstream  boundary  condition. 

For  potential  flow  F  =  0 . 

A  body  fitted  coordinate  system  is  introduced  by  the  following  general 
transformation: 

x  -  x(  C,n)  ,  r  =*  r(£,n)  ,  (4) 

where  the  coordinates  ( ?, n)  may  or  may  not  be  orthogonal.  Transforming  to 

(£,n)  coordinates,  Eq.  (1)  becomes 


where 


G  -  AG  -  2BG-  +  CBrc  +  DG,  =  -  E  , 

nn  n  £ 


A-  ^  [US)E-  , 


B  =  -  , 

Y  ’ 


«-■?  ■ 


D  -  ^  [Uc)s  -  <JB)„  +  , 
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2 

E  “  ^2~  F  , 

J  Y 

and 

J  =  Jacobian  =  (xr  -  x  r  1  * 

1  E  n  n  V 


do: 


(ii 


a 


2  2 

x  +  r 

n  n 


(12 


8  =  x_x  +  r_r 

5  n  5  n 


(13 


Y 


2  2 

x5  +  r5 


(14 


Body-Fitted  Coordinates 

To  produce  suitable  body-fitted  coordinates  an  analytically  defined 
C-grid  is  used,  as  illustrated  in  Fig.  1,  consisting  of  three  regions. 
The  transformations  are  as  follows: 

(I)  Orthogonal  wrap-around  grid,  0  <  £  <  £  . 


x  =  xfe(£)  -  n 

sin 

<KO 

(15 

r  *  rfe(£)  +  n 

cos 

♦U) 

(16 

where  x^,  r^ 

and 

as 

well  as  £  and 

n  are  defined  in  Fig.  2. 

(II)  Sheared 

grid, 

Si 

<  «  <  V 

dx  =  cos  <t>(£) d£ 

(17 

r  -  tb(C)  ♦ 

n  . 

(18 

We  note  that  for  Regions  I  and  II  to  be  compatible  the  junction  must  occur 
at  the  maximum  body  diameter. 
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(III)  Cartesian  grid,  ^ 

x  -  £  (19) 

r  =*  n  .  (20) 

The  metric  coefficients  x  ,  x  ,  etc.  are  obtained  analytically  from 

£  n 

Eqs.  (15)  through  (20). 


Potential  Flow  Boundary  Conditions 

For  the  potential  flow  case  Dirichlet  boundary  conditions  are  specified 
on  the  four  sides  of  the  computational  domain  consisting  of  Regions  I,  II  and 

III.  On  the  stagnation  line  (left  boundary)  \|>  =  0  and  hence 

G(0,n)  »  0  ,  0  <  n  <  Hy  .  (21) 

Since  the  body  and  wake  are  part  of  the  same  streamline  as  the  stagnation 
line,  the  condition  i|>  *  0  holds  which  leads  to 

f  -  j  r l  ,  0  <  5  <  ,  (22a) 


G(C,0) 


0 


£  <  £  <  £ 

T  R 


(22b) 


At  the  outflow  plane  (right  boundary)  the  flow  is  assumed  to  be  uniform  and 
parallel,  the  same  as  the  free  stream,  and  hence 

G(£r»t0  =  0  ,  0<ri<nu  •  (23) 

Finally,  the  outer  boundary  is  assumed  far  enough  removed  from  the  body  so 
that  free-streara  conditions  prevail.  Thus 

G ( £* r\j )  =  0  0  <  5  <  £r  •  (24) 


Frozen  Vorticity  Boundary  Conditions 

In  the  frozen  vorticity  case  the  computational  domain  consists  of  part 
of  Region  II  and  all  of  Region  III.  The  boundary  conditions  on  n  *  0  and 
n  ■  a  are  the  same  as  in  the  potential  flow  case.  On  the  left  or  inflow 
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boundary  the  stream  function  and  vorticity  are  specified  as  a  composite  of  the 
potential  flow  and  bound ary- layer  solutions,  viz. 

GUpv.h)  -  Gc(n) 

■  ,  0  <  n  <  .  (25) 

p(Svn)  -  Fc(n)J 

On  the  outflow  plane  (right  boundary)  parallel  flow  is  assumed,  as  in  the 
potential  flow  case,  which  leads  to 

cpUR,h)  =  0  ,  (26) 

plus  the  following  quadrature  relation  between  r  and  *{;: 

1  ^  2 

4>  =  j  f  udr  on  L  ,  0  <  n  c 


0 

where  from  Bernoulli's  equation  and  Eq.  (26) 


(27) 


u  =  /Cpo(«jO  * 


(28) 


The  numerical  details  of  evaluating  the  inflow  and  outflow  boundary  conditions 
will  be  given  later. 


Numerical  Algorithm 

The  main  features  of  the  numerical  method  used  to  solve  the  potential 
flow/frozen  vorticity  equation  are  as  follows: 

(1)  The  transformed  vorticity  equation,  Eq.  (5),  is  written  as  a  first- 
order  system. 

(2)  A  fourth-order  accurate  spline,  S^/i.O)  is  used  in  the  n  direction  to 
resolve  the  vortical  layer  with  as  few  nodal  points  as  possible. 

(3)  Second-order  accurate  finite  difference  formulas  are  used  in  the 


^-direction. 
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(4)  A  nonuniform  grid  in  %  and  n  is  generated  by  the  use  of  stretching 
functions . 

(5)  The  resulting  system  of  algebraic  equations  is  solved  by  SLOR  sweeping  in 
the  ^-direction. 

To  write  Eq.  (5)  as  a  first-order  system  we  define  the  following 
auxiliary  variable: 

H  -  Gn  .  (29) 

Then  Eq.  (5)  may  be  written 

Hn  -  AH  -  2BH?  +  CG^  +  DG^  -  -  E  .  (30) 

The  next  step  is  the  definition  of  the  following  spline  first  derivatives: 

*G  -  Gn  ,  (31) 

*H  -  Hn  .  (32) 

Then  the  governing  equations,  (29)  and  (30),  become 

lG  -  H  -  0  ,  (33) 

and 

JlH  -  AH  -  2BH?  +  CG??  +  DG?  +  E  »  0  .  (34) 

The  finite  difference  expressions  used  for  the  ^-derivatives  are  those  given 
by  Blottner  [7]  for  a  nonuniform  grid,  namely 
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*1,3  '  al,3°l,3  '  bl,3Hi,3  ‘  Ei.j 


where 


i*j  C, 


b.  .  *  A, 

i,J  i.J 


1 _  ( _ I _ +  _ i _ )C 

-  q-i  l?i+1  -  h  h~  w 


\  J  - 1 — I-? —  [(T^f- -  J +  +  Di  j)ci+i  j] 

H+i  H-l  H  H-l  1,:!  1  1,3  H+l  H  1,;1  1 

2B  . 

(H^.i  ~  W  '  Vl  •  <40> 

G  H 

The  unknowns  at  node  point  (i,j)  are  G,  l ,  H  and  l  .  Therefore  to 
complete  the  system  two  spline  relations  are  needed.  The  governing  equation 
has  been  written  in  first-order  form  so  that  the  same  spline  relation, 

S1 (4 , 0),  can  be  used  twice.  We  therefore  have  (i  subscript  understood). 

AAj$j_1  +  a2 A j _ j  +  BB j  j  +  ( 1  +  a)2*J  +  ccj*j+i  +  *j+i  *  0  .  (41) 

where  $  denotes  G  or  H,  and 

AA  =  x2q--y,-+-.-q\  ,  (42) 

j  An.—1  ( 1  +  a) 


2(1  -  o) ( 1  +  a) 
Anj-1° 

-  2(1  +  2 a) 

+  a)  a  ’ 


""j 


(45) 
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An.  -  tv  ,  -  n.  .  (46) 

J  3+1  j 

The  number  of  unknowns  at  node  point  (i,j)  is  reduced  from  four  to  two  by 

G  H 

solving  Eqs.  (33)  and  (37)  for  l  and  l.  then  substituting  into  the  two 
spline  relations.  The  following  two  tridiagonal  relations  are  thus  obtained: 

+  ”  +  +  O  +  o)  %  j 

•  <47) 


°2al,j-lGl,j-l  +  (AA3  +  +  <l  *  '’>\,JGl,j 

+  [BBj  4  U  +  A.jK.J  +  +  (CCj  +  bi,J+.)Hi,3« 


-  -  ‘Vj-l  -  (1  -  <’)2Rl,j  -  Ri,3+1 

By  defining  the  two-component  column  vector 


Eqs.  (47)  and  (48)  may  be  written  as  the  following  tridiagonal  matrix 
equation: 


Bi,jZi,j-l  +  Ai,jZi,j  +  Ci,jZi,j+l  "  Di,j  ’  2  <  j  <  N  » 


where  the  2x2  matrix  elements  A,  B  and  C  and  the  two-component  column 

A 

vector  D  are  obtained  from  Eqs.  (47)  and  (48). 

The  boundary  conditions  plus  a  two  point  spline  relation  at  each 
boundary  are  used  to  close  the  system  of  equations  on  line  i.  The  boundary 


conditions  are 
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Xl 


1  2 


2  rb. 


5  <  Cj; 
5  >  Ct 


(51) 


and 


Gi,N+l  "  °  * 


(52) 


The  two-point  spline  relation  used  here  is  Eq.  (16)  of  Ref.  8  which  in  the 


present  context  is 
An 


-  G 


An, 


H 


An.  H 

+  - -  £H 


An 


1  *H  -  0 


i  ,2  i,l  2  i,l  2  i  ,2  12  i,2  12  i,l 


(53) 


and 


% 


An. 


An2 

Gi,N+l  "  °i,N  "  ~  Hi,N  ~  ~T  Hi,N+l  +  *i,N+l  ”  XT'  *i,N  *  °  *  (54) 


Equations  (53)  and  (54)  are  fifth-order  accurate  in  An  and  were  found  to 

j 


yield  more  accurate  results  on  a  model  problem  [Eqs.  (47)  and  (48)  with 
constant  coefficients]  than  two-point  relations  of  lower  accuracy.  Upon 
elimination  of  l  using  Eq.  (34)  the  following  matrix  equations  result  at  the 
boundaries: 


Ai,lZi,l  +  Ci,lZi,2  “  Di,l  • 


(55) 


and 


Bi ,N+lZi,N  +  Ai,N+lZi,N+l  “  °i,N+l 


(56) 


The  singularity  at  the  tail  point  and  on  the  centerline  requires  special 

considerations.  As  can  be  seen,  Eq .  (5)  at  r  =  0  reduces  to  G  =0.  This 

n 

condition  replaces  spline  relation  (53)  and  is  written  as 
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A  special  form  of  Eq.  (48)  for  j  ■  2  is  also  necessary  because  b  is 

singular  at  n  ■  0.  Thus  Eq.  (37)  cannot  be  used  to  eliminate  £.  . .  Instead 

*  >  * 

the  two-point  spline  relation  Eq.  (53)  is  solved  for  Z^  ^  which  upon 
substitution  into  Eq.  (47)  yields, 


1  r  °i,i  ♦  (m2  -  +  {^r ♦ 1°2  *  (1  +  <’)2iai,2,Gi,2 


An 


*  (BB2  -  ^  +  (1  +  «)2]bi,2)Hl,2  *  ai,3Ci,3 


2  2 

+  (cc3  +  bi>3]Hi>3  -  -  [o  +  (1  +  a)  ]r1j2  -  Ri  3  .  (58) 

Along  any  line  i  =■  constant  the  set  of  block  tridiagonal  equations  (50) 
(55)  and  (56)  is  solved  by  L-U  decomposition  using  case  (i)  of  Ref.  9. 


Map  Junction  Lines 

At  map  junction  lines  Eq.  (37)  is  modified  to  account  for  the 
discontinuity  in  the  metric  coefficients  using  the  generalized  Chmielewski- 
Hoffman  method  of  Ref.  10.  The  author  has  found  that  ignoring  these 
discontinuities  leads  to  errors  as  large  as  37  percent  in  the  potential 
flow  pressure  coefficient. 

In  the  C-H  method  each  adjoining  domain  is  extended  one  step  into  the 
other  to  form  a  line  of  fictitious  nodes.  Then  the  vorticity  equation, 

Eq.  (5),  is  written  in  the  left  and  right  regions  at  the  junction.  Coupled 
with  a  condition  on  smoothness  of  derivatives  across  the  junction,  a  single 
spline-finite  difference  (SFD)  equation  can  be  derived  on  the  junction  line 
which  accounts  for  the  discontinuity  in  the  metric  coefficients. 
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The  numerical  strategy  used  here  to  derive  the  SFD  junction  equation  is 
as  follows: 

(1)  The  step  size  in  the  5-direction  is  assumed  to  be  constant  on  either  side 
of  the  junction  thus  simplifying  the  algebra. 


(2)  Central  differences  are  used  5-derivatives  of  G  using  the  fictitious 
nodal  values. 

(3)  One-sided  differences  are  used  for  thus  avoiding  the  necessity  of 
fictitious  nodal  values  for  H. 

The  governing  equation  at  the  map  junction  for  the  left  and  right  regions 
may  be  written  as 


£H  -  A(£)H  -  2B^H^  +  C^G^  +  D^G^£)  -  - 


(59) 


£H  -  A(r)H  -  2B(r)HpC)  +  C^G^  +  D(r)Gpr)  -  -  E(r) 

4  55  5 


(60) 


where  superscripts  l  and  r  denote  the  left  and  right  regions  respectively. 
The  smoothness  condition  merely  requires  x-derivatives  on  either  side  of  the 
junction  to  be  continuous.  This  condition  translates  to 


(jyJUM£)  -  (JyJ(rMr)  -  [(JyJU)  -  (JyJ(r)]H  . 


n'  5  ~5  ““'5'  ^“'5 

Performing  the  discretization,  as  outlined  above,  leads  to  an  equation  of 

the  same  form  as  Eq.  (37)  but  with  the  following  coefficients: 

*U)rT„  lUMr)  A  „(r)r,„  i(r).  (£) 


(61) 


1,J  «2  ’ 


(62) 


-  [2B»{  - 


,(£). (r) 


(r) 


,(r)T(*) 


i.J 


t 
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--*?  t- 


C^^fjy  l^l/r^G  +  C^fjy  )^l/^G 

i.r  yn'i.l  l.i  i-l.i  i,r  Vi.J  i,j  i 

(jy  +  [jy 

Jyn  i|j  i.j  yt  i.j  i.j 


where 


,  BU?(jy  ,  .  -  B^r?  (jy  J^L^V  , 

_  i_  r  hi:  Vi.J  i.j  1-1, J _ *.|J.. /Vi,,!  LJ-JiL, 

AS  (Jy  +  (jy 

1  n;i,J  i.j  1  Vi,j  i,j 

_U)rT  dDAv)  ^  (r)  /  dr)  Al) 

._2  r/  \(&)  (r)  /•  \(r)  (£)-|  * 


C^)  +  i.  d(JI)A£  , 

i,j  2  i,j 

p(«  In(0A- 

Ci,j  2Di,jA5  * 


Calculation  of  Inflow  and  Outflow  Conditions 

For  the  frozen  vorticity  case  the  conditions  on  the  inflow  boundary  are 
all  important  in  determining  the  solution.  As  already  mentioned  these 
conditions  are  a  composite  of  the  potential  flow  and  boundary-layer  solutions 
at  the  initial  line.  This  composite  is  calculated  as  follows: 

(1)  The  total  head  in  the  boundary  layer,  c  ,  in  the  law-of-the-wall  region 

p0 

has  been  observed  to  have  a  nearly  constant  derivative.  Thus  at  the  beginning 

of  this  constant  derivative  region,  usually  around  y+  *  120,  the  total  head 

is  extrapolated  linearly  to  the  wall  yielding  (cp  )  . 

'0  W 

(2)  Knowing  (c  )  ,  a  corresponding  wall  slip  velocity  is  computed  from 

W 

Bernoulli's  equation.  The  x-component  of  this  slip  velocity  is 


I 
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"w  -  ~  «wl  V  -  =p//2  ■ 

u  w 

where  c  is  the  wall  static  pressure.  The  x-component  of  velocity  in  the 
PW 

extrapolated  region,  0  <  y  <  y^ ,  is  computed  from  a  quadratic  polynomial 

since  u,.  as  well  as  u  and  (u  )  are  known, 
w  exc  y  ext 

We  note  that  the  extrapolation  of  (c  )  ,  as  given  in  Step  (1),  will 

^0  W 

guarantee  a  realistic  wake  solution  since  (cPq)w  >  0  and  far  downstream  in 
the  wake  the  centerline  velocity  is  given  by 

u(x,0)  =*  /(cp  )  . 

0  W 

(3)  The  modified  boundary-layer  velocity  profile  (with  wall  slip)  and  total 

head  profile  are  then  merged  smoothly  with  their  potential  flow  counterparts. 

In  the  case  of  the  velocity  profile  the  merging  occurs  smoothly  and 

naturally.  In  the  case  of  the  total  head  profile  the  merging  point  is  taken 

where  c  in  the  boundary  layer  just  becomes  unity.  The  total  head  in  the 
p0 

boundary  layer  as  y  increases  will  exceed  unity  because  the  y-component  of 
velocity  continues  to  increase. 

(4)  With  u  known  on  the  inflow  boundary,  denoted  by  u^y,  the  stream  function 
is  determined  by  numerical  integration  of 


The  integration  is  carried  out  using  the  trapezoidal  rule  formula. 


(4V>. 


(UW  +  UTO  )  * 
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where 


X  =  r^ 

AX  =  X  -  X 


j-1 


(69) 

(70) 


The  perburbation  stream  function  Gm.  is  then  calculated  from 

rv 


GFV  ^FV.  "  2  xj 

j  3  J 


(71) 


(5)  With  (c_  )  known,  the  reduced  vorticity  distribution,  Fw..(r),  is 
hq  FV  * v 


determined  from 


FV 


1 

2ru 


FV 


3r 


(72) 


where  the  derivative  of  fc_  1  with  respect  to  r  is  calculated  using  a 

p0  FV 


three-point  unequal  spacing  finite  difference  formula. 

The  procedure  described  above  is  preferable  to  that  of  Ref .  6  because 
fewer  approximations  are  made.  Here  the  only  altering  of  the  profile  is  in 
the  "slip  layer"  for  y+  <  120. 


The  outflow  stream  function  distribution  is  also  calculated  using 

R 


the  trapezoidal  rule,  viz. 


(aar) 


V 


j-i 


(V  .  +  urJ  » 

3-1  3 


(73) 


where  u  depend  on  iji  through  Eq.  (28).  Thus  for  each  integration  step 
R  R 


iteration  is  required  to  determine  iJ'r  during  which  u^  is  allowed  to  lag  one 

j  3 

cycle.  The  foregoing  procedure  is  much  simpler  and  produces  results  almost 


as  accurate  as  solving  for  >1/  from  a  two-point  boundary  value  problem  using 

R 


spline  discretization 
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II.  Results 


Computer  Codes 


The  main  advantage  of  the  frozen  vorticity  approach  to  the  axisymmetric 
strong  interaction  problem  over  its  competitors  is  the  high  degree  of 
automation  possible  thus  making  its  use  relatively  painless.  Two  computer 
programs  are  involved  in  the  calculation  and  the  key  to  the  automation  is 
the  communication  of  these  codes  with  each  other  through  automatically 
generated  data  files.  The  two  codes  are: 

AXFL02  —  Performs  potential  flow/frozen  vorticity  calculations 
using  the  method  given  in  Section  II  of  this  report. 

BL20  —  Performs  boundary-layer  calculations  using  the  Keller 
Box  method  in  conjunction  with  the  improved  algebraic 
turbulence  model  for  axisymmetric  bodies  given  in 


Ref.  11. 


The  steps  in  a  frozen  vorticity  strong  interaction  calculation  are  as 
follows : 

(1)  Using  AXFL02  in  the  potential  flow  mode  a  body  pressure  distribution  is 
generated.  Two  data  files  are  created,  one  for  BL20  containing  the  pressure 
distribution  and  body  curve  fit,  and  the  other  for  AXFL02  (frozen  vorticity 
mode)  containing  the  potential  flow  profile  on  the  initial  line. 

(2)  BL20  is  run  next  to  obtain  a  boundary-layer  solution  using  the  potential 
flow  pressure  distribution  from  AXFL02.  A  data  file  is  created  for  AXFL02 
which  contains  the  boundary-layer  profile  on  the  initial  line.  To  obtain 


profiles  of  u  and  c 

P0 


in  cylindrical  coordinates  a  double  interpolation 


procedure  is  performed. 
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(3)  AXFL02  is  now  run  in  the  frozen  vorticity  mode  to  compute  the  modified 
body  pressure  distribution  in  the  strong  interaction  region.  The  program  uses 
the  two  initial  profile  data  files  created  previously  and  forms  a  composite 
profile  following  the  procedure  described  in  Section  II. 

At  the  present  time  AXFL02  can  handle  only  two  mapped  regions  on  the 
body,  as  shown  in  Fig.  1.  Thus  bodies  with  cylindrical  mid  sections  and  flat 
noses  which  require  more  than  two  regions  for  proper  definition  are  excluded. 

Grid  Point  Distributions 

A  nonuniform  point  distribution  in  £  and  n  is  generated  using  one¬ 
dimensional  stretching  functions.  In  the  £ -direction,  a  two-sided  Vinokur 
distribution  [12]  is  used  on  each  body  segment  whereas  a  geometric 
progression  is  used  in  the  wake.  The  £  distribution  must  satisfy  the 
requirement  that  A?  on  either  side  of  the  map  junctions  is  the  same.  In 
the  n-direction  a  one-sided  Vinokur  stretching  function  is  used.  The 
stretching  functions  of  Vinokur  were  chosen  because  they  produce  a  grid 
with  a  uniform  truncation  error  independent  of  the  governing  differential 
equation  or  difference  algorithm. 

Although  the  same  £  distribution  is  used  in  the  potential  and  frozen 
vorticity  calculations,  the  n  distributions  differ  considerably.  In  the 
frozen  vorticity  case  for  proper  resolution  about  a  third  of  the  grid 
points  must  be  placed  in  the  thin  vortical  layer  which  is  the  same  thickness 
as  the  boundary  layer.  To  do  this  requires  a  rapidly  expanding  grid  in  n 
since  typically  ~  2  which  is  the  main  reason  splines  were  chosen  to 
approximate  flow  derivatives  in  n  rather  than  finite  differences. 
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Body  Curve  Fit 

An  accurate  body  curve  fit,  a  prerequisite  in  obtaining  an  accurate 
potential,  boundary  layer  and  frozen  vorticity  solution,  is  achieved  as 
follows:  With  the  body  shape  r^(x)  given  either  analytically  or  discretely, 
the  corresponding  arc  length  distribution  £  is  computed  using  the  spline 
formula 


where  the  two  cases  that  must  be  considered  are  shown  in  the  following  table 


case 

=  x(r  ) 

D 

“  rbU) 


range 


r;  > 1 

0  <  r  *  <  1 
b 


Table  1. 


u 

m 

h 

,  .2^/2 
(1  +  x’Z) 

x’x'  ' 

u 

Arb 

+  r'b> 

r'r" 
b  b 

u 

Ax 

Body  Arc  Length  Parameters. 


The  first  case,  x  =  xCr^),  is  appropriate  in  the  nose  region  of  a  blunt-nose 

body.  The  quantities  x' ,  x' ' ,  r'  and  r’’  required  in  the  £  calculation  are 

b  b 

computed  using  a  cubic  parametric  slope  spline. 

Since  m^  depends  on  £^,  Eq.  (72)  must  be  solved  at  each  step  by 
iteration.  Convergence  is  very  rapid  usually  requiring  about  four  cycles. 
The  use  of  Eq.  (72)  will  produce  a  £  distribution  accurate  to  about  four 
decimal  places  thus  assuring  that  the  derivatives  required  in  the  mapping 
are  accurate  and  smooth.  Sample  calculations  have  shown  that  Eq.  (72) 
produces  a  £  distribution  three  to  four  times  more  accurate  than  the  chord 
formula  in  regions  where  r'  is  large  and  changing  rapidly.  Once  £  is  known 


-24- 


14  September  1984 
GHH : lhz 


for  the  input  distribution,  then  the  various  derivatives  of  r^  for  the 
desired  £  distribution  are  determined  by  parametric  spline  interpolation. 

In  determining  dug/dC  for  the  boundary-layer  calculation,  an  uneven  three- 
point  difference  formula  is  used  because  it  has  been  found  to  produce  fewer 
wiggles  than  a  cubic  spline. 

Numerical  Solutions 

The  same  two  bodies  are  used  as  test  cases  in  this  report  as  were  used 
in  Ref.  6,  namely,  the  F-57  low-drag  body  of  Parsons  and  Goodson  [13]  and  the 
modified  spheroid  of  Patel  et  al.  [14].  Their  choice  was  dictated  by  the 
high  quality  of  experimental  data  available  for  each. 

The  grid  parameters  used  in  the  SFD  solutions  are  given  for  the 
potential  flow  and  frozen  vorticity  cases  in  Tables  2  and  3. 


F-57  and  Spheroid 

XR 

2.1 

*11 

2.0 

N5 

55 

k) 

35 

Z 

N 

21 

r\ 

An 

0.01119 

1 

Table  2.  Potential  Flow  Solution  Grid  Parameters. 
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F-57 

F-57 

Spheroid 

0.756 

0.707 

0.785 

2.1 

2.1 

2.1 

2.0 

2.0 

2.0 

28 

29 

28 

8 

9 

8 

41 

41 

41 

0.00055 

0.00055 

0.00055 

Table  3.  Frozen  Vorticity  Solution  Grid  Parameters. 


A  relaxation  factor  of  1.0  was  used  on  G  throughout  whereas  the 
relaxation  factor  on  lG  ranged  between  0.7  and  0.8.  Under-relaxation  on 
*  was  found  necessary  to  prevenc  divergence  of  the  calculation. 
Convergence  was  considered  accomplished  when  the  absolute  maximum  change 
in  H  was  less  than  10-^  which  in  both  potential  flow  and  frozen  vorticity 
runs  occurred  between  111  and  127  iterations  (note  that  the  number  of 
grid  points  in  each  case  was  nearly  the  same,  approximately  1150). 

Typical  CPU  times  on  a  VAX  11/782  computer  using  double  precision 
arithmetic  were  between  two  and  three  minutes. 

Boundary-layer  solutions,  as  already  mentioned,  used  the  same  £ 
spacing  as  the  potential  flow  calculations  plus  from  30  to  50  points 
nonuniformly  spaced  in  the  normal  direction.  Typical  CPU  times  were 


about  45  seconds . 
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To  determine  how  accurate  the  SFD  potential  flow  solutions  were  on  a 
55  x  21  grid,  the  Douglas-Neumann  program  with  80-90  surface  points  was  used 
as  the  standard.  Comparisons  of  body  pressure  distributions  for  the  F-57 
and  modified  spheroid  are  shown  in  Figs.  3  and  4.  Except  for  a  slight 
difference  near  the  nose  of  the  F-57.  the  results  are  indistinguishable. 

A  typical  composite  initial  velocity  profile  for  the  frozen  vorticity 

calculation  is  shown  in  Fig.  5,  in  this  case  for  the  F-57  at  xpy  =  0.756. 

The  corresponding  composite  total  head  profile  is  shown  in  Fig.  6  and  clearly 

exhibits  the  nearly  linear  region,  the  beginning  of  which  is  extrapolated 

to  the  wall.  The  slight  overshoot  of  c  at  the  edge  of  the  boundary  layer 

P0 

is  also  shown. 

The  one  disadvantage  of  the  frozen  vorticity  approach  is  that  the  user 
must  have  some  idea  where  to  begin  the  calculation,  i.e.,  where  to  choose 
Xpv-  Experience  has  shown  that  at  Reynolds  numbers  of  approximately  10^ 

(the  range  for  the  F-57  and  modified  spheroid)  that  the  strong  interaction 
region  begins  at  about  75  percent  of  the  chord  back  from  the  nose.  Thus 
the  question  naturally  arises  as  to  how  dependent  is  the  solution  on  the 
choice  of  xpy.  To  partially  answer  this  question  two  locations  were  chosen 
for  xpY  (0.707  and  0.756)  for  the  F-57  solutions.  Figure  7  gives  the 
comparison  of  body  pressure  coefficient  for  these  two  values  of  x^y.  The 
maximum  difference  in  Cp  between  the  two  solutions  is  0.0009  at  xg  =  0.81, 
or  about  6  percent. 

One  effect  evident  in  Fig.  7  is  the  upswing  in  cp  at  the  tail  which  is 
caused  by  the  numerical  scheme  attempting  to  simulate  a  stagnation  point. 

This  effect  is  of  course  absent  in  the  displacement  body  method  because  the 
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displacement  body  does  not  have  a  rear  stagnation  point.  To  eliminate  the 

rear  stagnation  point  in  the  frozen  vorticity  solution  a  cusped  tail  fairing 

was  arbitrarily  added  to  the  last  step  on  the  body  and  the  first  step  in  the 

wake.  This  fairing  was  taken  to  be  a  fourth-order  polynomial  for  which  r*' 

b 

was  continuous.  The  total  length  of  the  fairing  was  adjusted  to  force  r^' 
to  zero  at  the  junction  with  the  wake. 

Using  the  tail  fairing  approximation  frozen  vorticity  solutions  were 
obtained  for  the  two  test  cases.  The  results  for  body  pressure  distribution 
are  presented  in  Fig.  10  for  the  F-57  (xpy  =  0.756)  and  in  Fig.  11  for  the 
modified  spheroid  (xpy  =  0.785).  Comparisons  with  published  experimental 
data  are  given  in  each  figure.  For  the  F-57  overall  agreement  is  excellent. 
The  modified  spheroid  also  agrees  well  with  experiment  up  to  x  of  about 
0.95.  Beyond  this  point  the  prediction  is  too  high  with  a  maximum  error 
occurring  at  the  pressure  peak  (x  «  0.98)  of  some  25  percent.  This  dis¬ 
agreement  is  probably  attributable  to  the  22.4  degree  tail  angle  which  would 
make  the  tail  fairing  fairly  critical.  In  the  F-57  case,  the  tail  angle  is 
only  5.1  degrees,  almost  a  cusp. 


I 

« 


-28- 


14  September  1984 
GHH: lhz 


IV  -  Concluding  Remarks 

For  the  two  test  cases  examined  in  this  report  the  present  axisymmetric 
frozen  vorticity  procedure  has  been  found  to  give  pressure  distributions  in 
the  strong  interaction  region  of  acceptable  accuracy.  The  method  has  the 
advantage  over  its  displacement  body  rival  that  it  can  be  highly  automated 
and  requires  only  three  passes  to  obtain  the  final  results  —  a  potential 
flow  solution,  a  boundary-layer  solution  and  a  frozen  vorticity  solution. 

On  the  other  hand,  the  frozen  vorticity  method  suffers  from  the  disadvantage 
of  requiring  prior  knowledge  of  where  the  strong  interaction  region  begins. 

It  also  requires  a  cusped  fairing  to  be  placed  on  the  tail  of  the  body  to 
prevent  the  occurrence  of  a  rear  stagnation  point.  With  regard  to 
sensitivity  to  the  location  of  the  frozen  vorticity  initial  line,  preliminary 
results  show  an  error  of  six  percent  in  cp  with  a  five  percent  shift  in 
initial  line  location.  The  code  as  it  presently  exists  is  restricted  to 
smooth  bodies  which  can  be  described  by  two  segments. 
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Figure  7.  Effect  on  Pressure  Distribution  of  Initial  Plane  Location, 
F-57  Body. 


Figure  8.  Effect  on  Reduced  Vorticity  of  Initial  Plane  Locati 
F-57  Body. 
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Figure  9.  Effect  on  Outflow  Plane  Velocity  Profile  of  Initial  Plane 
Location,  F-57  Body. 
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